Noro-Prenkel scaling in short-range square well: A Potential Energy Landscape study 
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We study the statistical properties of the potential energy landscape of a system of particles 
interacting via a very short-range square- well potential (of depth — «o), as a function of the range of 
attraction A to provide thermodynamic insights of the Noro and Frenkel [ M.G. Noro and D. Frenkel, 
J.Chem.Phys. 113, 2941 (2000)] scaling. We exactly evaluate the basin free energy and show that 
it can be separated into a vibrational (A-dependent) and a floppy (A-independent) component. 
We also show that the partition function is a function of Ae'^"°, explaining the equivalence of the 
thermodynamics for systems characterized by the same second virial coefficient . An outcome of our 
approach is the possibility of counting the number of floppy modes (and their entropy). 
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Colloid and protein systems are often characterized 
by interactions on a length scale significantly smaller 
than the particle size. When the interaction range is 
less than 10% of the repulsive diameter, both termo- 
dynamic and dynamic properties are drastically differ- 
ent from the standard liquid behavior 1, 2, 3, 0, [3- 
An important characteristic of short-range spherical po- 
tentials V{r) is the independence of thermodynamic 
properties by the specific V(r) shape. In particular, 
the reduced second virial coefficient B2 = 62/62^ = 
-27r / (1 - e-'^^('^)) r'^dr/B^'^ (where ^ is the hard- 
sphere B2) appears to be a proper scaling variable which 
combines temperature T and potential shape. Based on 
the analysis of several available theoretical and numerical 
data, Noro and Frenkel (NF) have proposed a gener- 
alized law of corresponding states which states that all 
short-ranged spherically symmetric attractive potentials 
are characterized by the same thermodynamics proper- 
ties if compared at the same reduced density and 
This universal behavior is captured by the extreme limit 
of vanishing interaction range, i.e. by the celebrated 
Baxter sticky sphere model |7| , which, despite its known 
pathologies [3, has been successfully used to interpret 
experimental results in colloidal and proteins systems. 
The available accurate estimates of the critical point loca- 
tion and two-phase coexistence of the Baxter model ^3 
are becoming the reference for the thermodynamic be- 
havior of all short-range spherical potentials. 
In this Letter we study the NF empirical scaling in 
the potential energy landscape (PEL) thermodynamic 
framework, focusing on particles interacting via square- 
well (SW) potential of different well-width A, to provide 
a deeper understanding of the generalized law of corre- 
sponding states. We first show that for this model the 
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NF-law holds in a large range of densities (not only close 
to the liquid-gas critical point) and correctly reproduce 
the Baxter behavior in the limit of vanishing A. Then we 
calculate, with no approximations , the statistical proper- 
ties of the PEL [TH ll^ 0^ and their dependence on the 
well width. This study allows us to separate the basin 
free entropy in two components, respectively reflecting 
the exploration of the bond volume (A — dependent) and 
the exploration of the available volume at fixed bond dis- 
tance (A — independent) . This second contribution is the 
only one existing in the Baxter model. 
To generate equilibrated configurations, we perform 
event-driven molecular dynamics (EDMD) simula- 
tions for the SW potential, defined by a potential 
well of depth uq = — 1 and well width A. We investigate 
values of A/d ranging from 10~^ to 10""^, where c? = 1 is 
the hard-core particle diameter, chosen as unity of length. 
Pairs of particles are considered bonded when their rela- 
tive distance is between d and A. Density is expressed 
as packing fraction (f) — ^pd^, where p = N/L^, N being 
the number of particles and L the simulation box size. 
According to the NF law, the thermodynamic of the sys- 
tem is controlled by B2 and (p, at least in the proximity 
of the critical point. In the SW case B2 — I ^ Ae"°/^ in 
the limit of short well-width. 

Fig-fflshows the potential energy per particle U/N (panel 

a) and the compressibility factor Z = PV/NksT (panel 

b) as a function of cj) for several A values, at B'2 — —0.40. 
For A < 0.05, all curves collapse on the same curve con- 
firming the validity of NF scaling also outside the critical 
region. Similar behavior is observed for other different 
B2 values. The value at which the scaling breaks, i.e. 
A ~ 0.05, is in close agreement with the value at which 
the Baxter solution ceases to be a valid approximation 
for the SW model To confirm that the SW poten- 
tial approaches the Baxter model in the limit of vanish- 
ing A (and as a confirmation of our ability to equilibrate 
configurations at small A values) we compare in Fig. ^ 
the recently calculated radial distribution function g(r) 
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FIG. 1: (Color online) Top: (a) Energy per particle vs. </!> for 
various A at B2 = —0.40. (b) same as (a) for compressibility 
factor Z = PV/NkBT. (c) Radial distribution function for 
4> = 0.164 and = —0.92. The continuous line is from 
simulation of the Baxter model [T^ . dots are for a SW with 
A = 10-^ 



of the Baxter model with the the g{r) calculated for 
a SW system with A'= 10"'' at = -0.92. The agree- 
ment between the two results is perfect. The shoulders 
and the (5-functions indicating the presence of clusters 
with well defined structure |la | are reproduced both in 
their location and their height. 

The fact that for A < 0.05 the state of the system is 
controlled by the values of cj) and of Ae'^"" , suggests that 

— at fixed — the configurational part of the parti- 
tion function Z must be, in leading order, a function of 
Ae^"o. In the PEL formalism 11, 12,13] the phase space 
is decomposed into basins of attractions of the local min- 
ima of the PEL — the so-called inherent structures (IS) 

— and the partition function is written as a sum over 
all landscape basins, each of them weighted by the local 
minimum Boltzmann factor and by the basin free-energy. 
In the SW case, we associate an IS to a specific bonding 
pattern 113 (thus, the IS energy is expressed by the total 
number of bonds Nb) and the basin free-energy to the 
logarithmic of the (multidimensional) volume Q which 
can be explored in configuration space without breaking 
or reforming bonds. In this respect, the conceptual oper- 
ation of thinking Z as a sum over all possible PEL basins 
is equivalent to expressing Z as a sum over all possible 
bonding pattern Tq . Since the Boltzmann factor weight- 
ing the probability of each configuration is e^""^"" , scaling 
in the variable Ae^"° implies that the number of states ft 
sampled by each fixed bonding pattern must scale with 
A as A^"" , i.e. each bond is independent and contributes 
a quantity of order A to il. This is the hypothesis that 
we test next. 
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FIG. 2: (Color online) Top: Mean square displacement of a 
particle harmonically bounded to a reference point and in- 
teracting via a SW with a particle located in the origin for 
different A values. The probability distribution of the particle 
is schematically shown above the plot for three representative 
values of A: particle completely delocalized (left), particle 
localized by the bond (center) and particle localized by the 
harmonic force (right). Ai and A2 are discussed in the text. 
Bottom: Mean square displacement 2/3A(^^(ri — fi)'^)\/N 
(labelled 2(3X{r^) for semplicity in the figure) for a SW sys- 
tem with (j) = 0.30 and = —0.69 vs. A , for different A 
values. The Einstein limit is shown as a dashed line. 



To precisely estimate we implement a generalization 
of the thermodynamic integration technique first used 
by Ladd and Frenkel to evaluate the free energy of the 
hard sphere crystals 0. In this technique, the system 
Hamiltonian Hq is complemented by an harmonic term 
A^j(ri — fi)^ centered around a reference configuration 
(fi, . . . jTn) of strength A acting on each particle. The 
essence of this method lies in the possibility of perform- 
ing a thermodynamic integration over A from a state of 
known free energy to the state of interest (A = 0). When 
A is very large, particles harmonically oscillate around the 
reference state and the system behaves like a collection 
of 3N harmonic oscillators, i.e. an Einstein solid. If the 
unperturbed system is identified with a specific bonding 
pattern and Hq is defined as zero if the configuration has 
the correct bonding pattern and infinity if the bonding 
pattern has changed (i.e. if a bond is broken or a new 
bond is formed), then the result of the thermodynamic 
integration provides an exact measure of the basin free 
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energy. Since the exploration of space within the fixed 
bond-pattern basin takes place on a flat energy surface 
the only contribution to the basin free energy is entropic. 
The basin entropy per particle in units of ks, CTbasin, can 
be formally written as 



CTba 



-/3/e(T,Aoo) + 



0')ArflnA (1) 



where /e (T, A) is the free energy of 3A^ harmonic oscil- 
lators coupled by an elastic constant A, {■)\ an ensemble 
average at a fixed value of A and Aoo is any value of A 
such that the harmonic contribution in dominant as com- 
pared to ifo and thus (A X]i(ri — fi)^)A ~ ^NksT. From 
a technical point of view we evaluate (A^j(ri — fi)^)^ 
for thirty or more A values between Aqo and via MC 
simulations rejecting all moves which modify the bond- 
pattern. 

To guide the interpretation of the numerical results it is 
helpful to examine the behavior of two particles bonded 
by a square-well in three dimensions. Fixing particle 1 at 
the origin (to neglect the trivial center of mass degrees of 
freedom) and selecting an Einstein reference site acting 
on particle 2 in an arbitrary position r2 = (x, 0, 0) (with 
d < X < d + A), the resulting ((r2 — ^2)^)x is 

A((r2-f2)2)A = -^logZA(/3), (2) 

where the generating (or partition) function Zx((3) is 



Zx{p) = 2ti 



d+A I3\(r+xf 



(e 



1) 



2/3A 



rdr (3) 



The resulting A dependence of 2/3A((r2 — r2)^)A is shown 
in Fig.|3t for three different values of A. At large A (the 
harmonic hmit), the quantity 2/3A((r2 — f2)^)A goes to 
three, the total number of degrees of freedom. On de- 
creasing A, the function shows a two-step decay to zero 
with an intermediate plateau at the value two. The two 
cross-overs (from 3 — > 2 and from 2-^0) are taking place 
at values of Ai « (2/A)^ and A2 ~ (2/d)^. To interpret 
this behavior we recall that for very large A (A > Ai) 
confinement is provided by the harmonic potential. For 
A2 < A < Ai, confinement of the harmonic potential has 
become larger than A and the bond distance becomes 
the relevant quantity controlling the mean square dis- 
placement. For even smaller A (A < A2) confinement is 
provided by the finite bond volume (a corona of width A 
and inner surface Aird'^). The two-step cross-over makes 
it possible to count the number of modes which are re- 
lated to exploration of the bond width and the number 
of modes which are related to exploration of space at 
fixed interparticle distance. Indeed, the first cross-over 
is A-dependent while the second one is A-independent. 
Hence, the FL method allows one not only to evalu- 
ate the total change in entropy but also to count (and 
separate !) the number of modes which are related to 
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FIG. 3: (Color online) Top: Number of state sampled Q}'^ = 
exp((Tbasin) VS. A for diflterent state points. The straight lines 
are single parameter fit of the form ylA""" where rib is the 
number of bonds per particle per particle, i.e. = N^/N. 
The inset shows the the location of the corresponding state 
points in the — phase diagram. It also show the perco- 
lation line (dashed line) and the liquid-gas coexistence from 
RefE Bottom: Symbols are 2l3\{Y^-{vi - fif)x/N (labelled 
2/3A(r^) for semplicity in the figure) vs. A for A = 10~^; The 
dashed horizontal line are 2> — rib. Note that the intermediate 
A plateau is well reproduced by the fraction of fioppy modes 

ff = S- Tib. 



the exploration of the bond-distance (vibrational modes) 
from the modes which are related to the exploration of 
the volume with rigid bonds (floppy modes). The be- 
haviour of 2/3A(X]i(ri — ri)^)\/N for a bond configura- 
tion of 200 particles (a typical equilibrium configuration 
produced by MD simulation) for different values of A 
is shown in Fig. The bonding network is identical 
for all A values. As in the two-particles example, the 
shape of the curve shows two parts, one A-independent 
and one A-dependent. As from Eq. ^ the area under 
the A(^j(ri — ViY)\/N vs In A curve is a measure of 
fbasinj the entropy associated to the exploration of space 
at fixed bonding pattern. The A-independent and the 
A-dependent parts of A(^j(ri — ri)^)A give rise to two 
different contribution to dbasin, which we can identify as 
the floppy j20| A-independent — the dark shaded region 
in the cartoon in Fig. |2K — [Nuno-ppy = InfJfloppy) and 



the vibrational A-dependent part — the light coloured 
region— (iVavib = Infivib)- 

According to the PEL picture of the NF law, — 
gO-baain niust be proportional to A"'', where rib = N^/N 
is the number of bonds per particle. This prediction 
can be put under severe test, since both A and rib are 
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a-priori known, by comparing the calculated values for 
gO-baain "v^jth A"*" . Tlils comparisou is reported in Fig.|3t 
for several typical bonding configurations with differ- 
ent bonding patterns, extracted from equilibrium sim- 
ulations at different (p and T, encompassing the range 
f.03 < 716 < 2.6f. In all studied cases, fibasin ~ 
with extreme accuracy. The validity of such relationship 
suggests that for short-range SW, each bond acts as an 
independent unit (so that the vibrational entropy of the 
bonding pattern coincides with the sum of all the bond 
entropies). In the PEL language, the NF scaling is an 
expression of the independence of the bonds. In this re- 
spect, proper scaling between different potential shapes 
is predicted when the vibrational free energy of a bonded 
pair is chosen as a proper scaling variable. The possibil- 



10-' ' lo" lo' 

fr(3-N,) 

FIG. 4: (Color online) Entropy of the floppy modes (crfloppy) 
vs fraction of floppy modes (// — 3 — Nt/N). The straight 
line is a power low with exponent 3, i.e. o-fjoppy oc /|. 

ity of separating in a precise way flvih and fifloppy (see 
the cartoon in Fig|2K), allows us to evaluate also the (A- 
independent) volume in configuration space sampled by 
the a specific bonding pattern when all bond distances 
are fixed. This volume corresponds to the free rolling mo- 
tions of the particles with no bond-breaking of -forming 
and it is essentially the basin volume accessible to the 



Baxter model. It is interesting to investigate the depen- 
dence of this quantity on the number of bonds, since one 
expect that on increasing the connectivity, the entropy 
of the floppy modes should decrease. For the short-range 
S W under consideration the fraction of floppy modes (the 
height of the plateau in 2/3A(^-(ri — fi)'^)\/N) is found 
to be = 3 — ni, (see Fig. ^p) , a value consistent with 
the existence of Nf, independent vibrational degrees of 
freedom. Hence, the total floppy entropy (the area under 
the curve - see Eq. should vanish close to the point 
when rif, w 3 (i.e. each particle is involved in average 
in six bonds). This expectation is indeed born out in 
the calculated //-dependence of crfloppy shown in Fig. 0] 
Unexpectedly, the floppy entropy goes to zero as a power 
law CTfloppy fx The value of the exponent suggests that 
not only the number of modes but also the floppy entropy 
per mode vanishes with a power-law close to f f — 0. It 
is interesting to observe that at the state point where 
Nb « 3A^, the full entropy of the system would be given 
only by the logarithm of the number of topologically dif- 
ferent bonding patterns. 

In conclusion, we have shown that the empirical NF-law 
of corresponding states can be explained in the rigorous 
PEL thermodynamic formalism as arising from the addi- 
tive contribution of the bonds to the basin entropy. An 
unexpected outcome of this study is a new methodol- 
ogy to separately evaluate the fraction of floppy modes 
and their entropy for any specific bonded configuration, 
a method which can be used in studies of the rigidity of 
hard particle systems [2ll |. 
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